1 Описание данных
В данной работе изпользуются данные о ценности жилья в городе Бостон в 1970-1980-х годах (пакет MASS).
Основной акцент будет сделан на то, как средняя стоимость домов (переменная medv, измеренная в тысячах долларов), занимаемых владельцами, зависит от различных факторов.
Авторы работы измерили значения 14 различных параметров (включая medv) для 506 домов.
Цель — построить линейную модель для предсказания стоимости домов.
Используемые обозначения:
- Median value of owner-occupied homes in $1000s (medv)
- Per capita crime rate by town (crim)
- Proportion of residential land zoned for lots over 25,000 sq.ft (zn)
- proportion of non-retail business acres per town (indus)
- Charles River dummy variable (= 1 if tract bounds river; 0 otherwise) (chas)
- Nitrogen oxides concentration (parts per 10 million) (nox)
- Average number of rooms per dwelling (rm)
- Proportion of owner-occupied units built prior to 1940 (age)
- Weighted mean of distances to five Boston employment centres (dis)
- Index of accessibility to radial highways (rad)
- Full-value property-tax rate per/$10,000 (tax)
- Pupil-teacher ratio by town (ptratio)
- 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town (black)
- Lower status of the population (percent) (lstat)
Примечание в конце анализа я осознал, что переменная rad, в целом, является фактором, но у него большое количество уровней, поэтому я решил не переделывать, подумав, что это не сильно повлияет на анализ.:)
Давайте посмотрим, как выглядят данные.
house_cost <- Boston
house_cost$chas <- factor(house_cost$chas)
house_cost$chas <- relevel(house_cost$chas, ref = '0')
str(house_cost) # Все ли нормально открылось## 'data.frame': 506 obs. of 14 variables:
## $ crim : num 0.00632 0.02731 0.02729 0.03237 0.06905 ...
## $ zn : num 18 0 0 0 0 0 12.5 12.5 12.5 12.5 ...
## $ indus : num 2.31 7.07 7.07 2.18 2.18 2.18 7.87 7.87 7.87 7.87 ...
## $ chas : Factor w/ 2 levels "0","1": 1 1 1 1 1 1 1 1 1 1 ...
## $ nox : num 0.538 0.469 0.469 0.458 0.458 0.458 0.524 0.524 0.524 0.524 ...
## $ rm : num 6.58 6.42 7.18 7 7.15 ...
## $ age : num 65.2 78.9 61.1 45.8 54.2 58.7 66.6 96.1 100 85.9 ...
## $ dis : num 4.09 4.97 4.97 6.06 6.06 ...
## $ rad : int 1 2 2 3 3 3 5 5 5 5 ...
## $ tax : num 296 242 242 222 222 222 311 311 311 311 ...
## $ ptratio: num 15.3 17.8 17.8 18.7 18.7 18.7 15.2 15.2 15.2 15.2 ...
## $ black : num 397 397 393 395 397 ...
## $ lstat : num 4.98 9.14 4.03 2.94 5.33 ...
## $ medv : num 24 21.6 34.7 33.4 36.2 28.7 22.9 27.1 16.5 18.9 ...plot(house_cost) # скатерплот чтобы посмотреть структуру данных
2 Построение полной линейной модели
Для начала проведем стандартизацию переменных. Также на первом этапе не будем учитывать взаимодействие предикторов.
house_cost_scale <- as.data.frame(sapply(house_cost[-4], scale))
house_cost_scale$chas <- house_cost$chas
house_cost_scale$chas <- relevel(house_cost_scale$chas, ref = '0')Теперь построим полную линейную модель
mod_full <- lm(medv ~ ., data = house_cost_scale)
summary(mod_full)##
## Call:
## lm(formula = medv ~ ., data = house_cost_scale)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.69559 -0.29680 -0.05633 0.19322 2.84864
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.020206 0.023835 -0.848 0.396976
## crim -0.101017 0.030737 -3.287 0.001087 **
## zn 0.117715 0.034811 3.382 0.000778 ***
## indus 0.015335 0.045871 0.334 0.738288
## nox -0.223848 0.048126 -4.651 4.25e-06 ***
## rm 0.291056 0.031928 9.116 < 2e-16 ***
## age 0.002119 0.040430 0.052 0.958229
## dis -0.337836 0.045666 -7.398 6.01e-13 ***
## rad 0.289749 0.062813 4.613 5.07e-06 ***
## tax -0.226032 0.068912 -3.280 0.001112 **
## ptratio -0.224271 0.030796 -7.283 1.31e-12 ***
## black 0.092432 0.026662 3.467 0.000573 ***
## lstat -0.407447 0.039378 -10.347 < 2e-16 ***
## chas1 0.292128 0.093679 3.118 0.001925 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.516 on 492 degrees of freedom
## Multiple R-squared: 0.7406, Adjusted R-squared: 0.7338
## F-statistic: 108.1 on 13 and 492 DF, p-value: < 2.2e-163 Диагностика модели
3.1 Проверка влиятельных наблюдений
Ни одно значение не превышает условного порога в 2 единицы. Влиятельных наблюдений нет (Но, если честно, то видно, что есть ряд наблюдений, которые могут смещать оценки коэффициентов модели, так как для них значения расстояния Кукa больше).
ggplot(mod_full_diag, aes(x = 1:nrow(mod_full_diag), y = .cooksd)) +
geom_bar(stat = "identity") +
geom_hline(yintercept = 2, color = "red") +
xlab('Номер наблюдения') +
ylab('Расстояние Кука')
3.2 Проверка линейности взаимосвязи и гетероскедастичности
График распределения остатков выглядит не очень хорошо. Достаточно большое количество наблюдений находятся за пределами двух стандартных отклонений, а также отчетливо виден паттерн в остатках. Все это свидетельствует о наличии нелинейности во взаимосвязи, а также о гетероскедастичности (непостоянство дисперсии).
gg_resid <- ggplot(data = mod_full_diag, aes(x = .fitted, y = .stdresid)) +
geom_point() +
geom_hline(yintercept = 0) +
geom_smooth(method = "lm") +
geom_hline(yintercept = 2, color = "red") +
geom_hline(yintercept = -2, color = "red") +
xlab('Предсказанное значение') +
ylab('Oст. модели')
gg_resid
3.3 Проверка на нормальность распределения предсказанных значений
Квантильный график выглядит не очень хорошо, сказать, что стандартизованные остатки распределены нормально, нельзя.
qqPlot(mod_full_diag$.stdresid, xlab = 'Квантили нормального распределения', ylab = 'Квантили распределения остатков модели')
## [1] 369 3723.4 Проверка отсутствия коллинеарности предикторов
Если обратиться к структуре датасета, то можно заметить наличие корреляций между предикторами. Поэтому в данном случае проверка на мулькиколлинеарность особенно актуальна.
Один из способов проверки модели на мультиколлинеарность предикторов это использование VIF (variance inflation factor).
Если предиктор имеет значение VIF выше 2, то его следует исключать из модели.
Если таких предикторов несколько, то прменяется пошаговый алгоритм: расчет VIF, удаление предиктора с максимальным VIF, расчет VIF для обновленной модели до тех пор, пока все значения не будут меньше порога.
vif(mod_full)## crim zn indus nox rm age dis rad
## 1.792192 2.298758 3.991596 4.393720 1.933744 3.100826 3.955945 7.484496
## tax ptratio black lstat chas
## 9.008554 1.799084 1.348521 2.941491 1.073995В даном случае легко можем видеть наличие мультиколлинеарности в нашей модели. В дальнейшем ряд предикторов надо будет убрать из модели.
4 Предсказания модели
Построим график предсказаний стоимости от переменной, которая обладает наибольшим по модулю коэффициентом. В данном случае это переменная lstat.
Мы не можем одновременно учесть изменчивость всех предикторов. Поэтому как правило выбирается один, который интересует вас больше всего (в нашем случае lstat). И относительно него создается тестовый датасет, где целевой предиктор принимает значения от минимального до максимального, а все остальные предикторы представлены своими средними значениями.
test_data <- data.frame(
lstat <- seq(min(house_cost_scale$lstat), max(house_cost_scale$lstat), length.out = 90),
crim <- rep(mean(house_cost_scale$crim), 90),
zn <- rep(mean(house_cost_scale$zn), 90),
indus <- rep(mean(house_cost_scale$indus), 90),
nox <- rep(mean(house_cost_scale$nox), 90),
rm <- rep(mean(house_cost_scale$rm), 90),
age <- rep(mean(house_cost_scale$age), 90),
dis <- rep(mean(house_cost_scale$dis), 90),
rad <- rep(mean(house_cost_scale$rad), 90),
tax <- rep(mean(house_cost_scale$tax), 90),
ptratio <- rep(mean(house_cost_scale$ptratio), 90),
black <- rep(mean(house_cost_scale$black), 90),
chas <- rep('0', 90)
)
Predictions <- predict(mod_full, newdata = test_data, interval = 'confidence')
MyData <- data.frame(test_data, Predictions)
Pl_predict <- ggplot(MyData, aes(x = lstat, y = fit)) +
geom_ribbon(alpha = 0.2, aes(ymin = lwr, ymax = upr)) +
geom_line() +
ggtitle("Множественная модель") +
xlab('Низкий статус населения, sd') +
ylab('Предсказанное значение цены за дом, sd')
Pl_predict 
Полученная модель далека от идеальной, как было показано выше, также мы не убрали предикторы, которые незначимо вдияют на зависимую переменную. Она уже может делать, какие-то предсказания, но ниже мы попытаемся сделать ее лучше.
5 Дополнительная часть
Вернемся к нашей полной модели и вспомним результаты проверки на коллинеарность.
vif(mod_full)## crim zn indus nox rm age dis rad
## 1.792192 2.298758 3.991596 4.393720 1.933744 3.100826 3.955945 7.484496
## tax ptratio black lstat chas
## 9.008554 1.799084 1.348521 2.941491 1.073995Будем пошагово удалять из модели по одному предиктору с наибольшим vif до тех пор, пока все значения не будут меньше порога (2).
mod2 <- update(mod_full, .~. - tax)
vif(mod2)## crim zn indus nox rm age dis rad
## 1.791940 2.184240 3.226015 4.369271 1.923075 3.098044 3.954446 2.837494
## ptratio black lstat chas
## 1.788839 1.347564 2.940800 1.058220mod3 <- update(mod2, .~. - nox)
vif(mod3)## crim zn indus rm age dis rad ptratio
## 1.785343 2.183394 2.872809 1.904013 2.875130 3.641492 2.533616 1.598944
## black lstat chas
## 1.339554 2.927273 1.057571mod4 <- update(mod3, .~. - dis)
vif(mod4)## crim zn indus rm age rad ptratio black
## 1.765881 1.758636 2.517520 1.879925 2.423551 2.507024 1.530992 1.339553
## lstat chas
## 2.926111 1.056840mod5 <- update(mod4, .~. - lstat)
vif(mod5)## crim zn indus rm age rad ptratio black
## 1.727110 1.751140 2.489395 1.307312 2.038976 2.497242 1.529761 1.304899
## chas
## 1.053402mod6 <- update(mod5, .~. - nox)
vif(mod6)## crim zn indus rm age rad ptratio black
## 1.727110 1.751140 2.489395 1.307312 2.038976 2.497242 1.529761 1.304899
## chas
## 1.053402mod7 <- update(mod6, .~. - rad)
vif(mod7)## crim zn indus rm age ptratio black chas
## 1.372553 1.738620 2.229252 1.281667 2.026045 1.369726 1.251785 1.052519mod_good <- update(mod7, .~. - indus)
vif(mod_good)## crim zn rm age ptratio black chas
## 1.342583 1.680416 1.216701 1.668190 1.345528 1.209887 1.043732Итоговая модель этого шага (в нестандартизованном виде):
- При значении факторной переменной chas = 0
medv = -7.06 - 0.11 * crim - 0.003 * zn + 7.06 * rm - 0.04 * age - 0.93 * ptratio + 0.015 * black
- При значении факторной переменной chas = 1
medv = -3.54 - 0.11 * crim - 0.003 * zn + 7.06 * rm - 0.04 * age - 0.93 * ptratio + 0.015 * black
5.1 Этап 1. Поиск оптимальной модели
Мы можем оставить модель в таком виде, а можем попробовать оставить только те предикторы, которые значимо влияют на стоимость домов (большое число предикторов обычно приводит к оверфитингу модели).
5.1.1 Пошаговый отбор предикторов по значимости
В этой работе я буду использовать алгоритм отбора предикторов backward selection (он же backward elimination). В качестве критерия отбора буду использовать частный F-тест.
drop1(mod_good_unstd, test = 'F')## Single term deletions
##
## Model:
## medv ~ crim + zn + chas + rm + age + ptratio + black
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 15509 1747.9
## crim 1 323.2 15832 1756.3 10.3772 0.0013592 **
## zn 1 1.6 15511 1745.9 0.0506 0.8220598
## chas 1 386.8 15896 1758.3 12.4206 0.0004638 ***
## rm 1 10225.3 25735 2002.1 328.3350 < 2.2e-16 ***
## age 1 411.9 15921 1759.1 13.2258 0.0003049 ***
## ptratio 1 1514.7 17024 1793.0 48.6365 9.843e-12 ***
## black 1 781.1 16290 1770.7 25.0805 7.646e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod_good_unstd <- update(mod_good_unstd, .~. - zn)
drop1(mod_good_unstd, test = 'F')## Single term deletions
##
## Model:
## medv ~ crim + chas + rm + age + ptratio + black
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 15511 1745.9
## crim 1 329.2 15840 1754.5 10.591 0.0012133 **
## chas 1 390.0 15901 1756.5 12.546 0.0004342 ***
## rm 1 10413.2 25924 2003.8 335.006 < 2.2e-16 ***
## age 1 514.5 16025 1760.4 16.553 5.498e-05 ***
## ptratio 1 1605.3 17116 1793.7 51.644 2.440e-12 ***
## black 1 780.3 16291 1768.8 25.103 7.559e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Таким образом, мы оставили только предикторы со значимым влиянием на зависимую переменную.
5.1.2 Диагностика модели
Мы уже смотрели на ряд графиков для диагностики нашей модели, но помимо них для множественных моделей необходимо строить графики от предикторов, не вошедших в модель. В нашем случае видно, что в модели, в общем-то, нет неучтенных зависимостей, соответственно, можно не возвращать предикторы, убранные из модели на предыдущих этапах анализа
res1 <- gg_resid + aes(x = zn)
res2 <- gg_resid + aes(x = indus)
res3 <- gg_resid + aes(x = nox)
res4 <- gg_resid + aes(x = dis)
res5 <- gg_resid + aes(x = rad)
res6 <- gg_resid + aes(x = tax)
res7 <- gg_resid + aes(x = lstat)
grid.arrange(res1, res2, res3, res4, res5, res6, res7, nrow = 4)
mod_good_unstd_diag <- data.frame(fortify(mod_good_unstd), house_cost[, c(1, 4, 6, 7, 11, 12)])
ggplot(mod_good_unstd_diag, aes(x = 1:nrow(mod_good_unstd_diag), y = .cooksd)) +
geom_bar(stat = "identity") +
geom_hline(yintercept = 2, color = "red") +
xlab('Номер наблюдения') +
ylab('Расстояние Кука')
gg_resid <- ggplot(data = mod_good_unstd_diag, aes(x = .fitted, y = .stdresid)) +
geom_point() +
geom_smooth(method = "lm") +
geom_hline(yintercept = 0) +
geom_hline(yintercept = 2, color = "red") +
geom_hline(yintercept = -2, color = "red") +
xlab('Предсказанное значение') +
ylab('Стандартизованные остатки модели')
gg_resid
qqPlot(mod_good_unstd_diag$.stdresid, xlab = 'Квантили нормального распределения', ylab = 'Квантили распределения остатков модели')
## [1] 369 372Анализ показал, что, если неприятности и исправились, то не очень сильно. Модель вышла не идеальная. Может быть, в данном случае помогла бы какая-то трансформация переменных (зачастую не желательна, так как дальше сложно интерпретировать результаты) или применение других статистических методов, о которых, вероятно, речь пойдет в будущих проектах. Также попытки учесть какие-то взаимодействия переменных внутри итоговой модели не привели к успеху, наоборот график расстояний Кука становился только хуже (не привожу здесь). В первом приближении мы можем попытаться предсказать, какие аспекты надо улучшить, чтобы максимизировать цену за дом.
5.2 Этап 2. Рекомендации заказчику
Итоговая модель:
- При значении факторной переменной chas = 0
medv = -7.06 - 0.11 * crim + 7.05 * rm - 0.04 * age - 0.92 * ptratio + 0.015 * black
- При значении факторной переменной chas = 1
medv = -3.54 - 0.11 * crim + 7.05 * rm - 0.04 * age - 0.92 * ptratio + 0.015 * black
Мы можем видеть 6 параметров, которой так или иначе влияют на цену домов. Вот несколько замечаний, которые, на мой взгляд, должны максимизировать цену за продаваемый дом:
- Участки, граничащие, с рекой в среднем стоят на 3.52 тысячи долларов больше
- При снижении уровня преступности на душу населения на одну единицу цена за продаваемый дом возрастает на 0.11 тысяч долларов
- В среднем каждая дополнительная комната увеличивает цену за дом на 7.05 тысяч долларов
- При уменьшении доли жилых домов, построенных до 1940 г. на одну единицу цена за дом увеличивается на 0.04 тысячи долларов. Этот параметр можно было бы вовсе не учитывать, если строить полностью новый район
- При уменьшении соотношения учеников и учителей в городе цена за дом возрастает на 0.92 тысячи долларов
- Параметр black отражает долю чернокожего населения в городе, но на него вряд ли можно будет повлиять, так что рекомендаций по этому поводу не будет
Таким образом в идеальном районе:
- Участки граничат с рекой
- Низкий уровень преступности
- Большое количество комнат в домах (в данном случае до 8, так как в данных нет более высокого показателя)
- Было бы хорошо, если бы количество учеников в классе было минимально, например, три
Примерная цена за дом при выполнении вышеуказанных условий может составить 50.1 тысяч долларов.